CASE STUDY 3: SEARCH FOR THE UNUSUAL CLUSTER IN THE PALINDROMES

Question

In this paper, we will search for unusual clusters of complementary palindromes. The overarching research question is: “How do we find clusters of palindromes? How do we determine whether a cluster is just a chance occurrence or a potential replication site? Based on our analysis, we will then provide recommendations to biologists who are about to start experimentally searching for the origin of replication.

Setup

locations <- read.table('hcmv-25kgjn1-1rfrtkc.txt', header=TRUE)$location  # Original
health <- read.csv('RAW_DATA-2iwcznn-2kr2xw0.csv', header=TRUE)  # Additional
N <- 229354  # Base pairs
n <- 296  # Palindromes

Scenario 1: Random Scatter

To begin, pursue the point of view that structure in the data is indicated by departures from a uniform scatter of palindromes across the DNA.

Of course, a random uniform scatter does that mean that palindromes will be equally spaced as milestones on a freeway. There will be some gaps on the DNA where no palindromes occur, and there will be some clumping together of palindromes.

To look for structure examine the locations of the palindromes, the spacing between palindromes, and the counts of palindromes in non overlapping regions of the DNA. One starting place might be to see first how random scatter looks by using a computer to simulate it.

A computer can simulate 296 palindrome sites chosen at random along a DNA sequence of 229,354 bases using a pseudo random number generator. When this is done several times, by making seller sets of simulated palindrome locations, then the real data can be compared to the simulated data.

set.seed(0)
color <- 'red'
# Generate 3 samples from the uniform distribution with the same size and bounds as our data
samples=list(sort(runif(n, min=0, max=N)), sort(runif(n, min=0, max=N)), sort(runif(n, min=0, max=N)))
# Dot plot of locations of palindromes in original data and uniform scatter
title1 <- 'Locations of Palindromes'
title2 <- c(title1,'(Simulated)')
x.axis <- 'Base Pair'
symbol <- 3
stripchart(locations, pch=symbol, col=color, main=title1, xlab=x.axis)

for (sample in samples) {
  stripchart(sample, pch=symbol, main=title2, xlab=x.axis)
}

# Additional dot plot of locations of palindromes in original data and uniform scatter
dotchart(locations, color=color, main=title1, xlab=x.axis)

for (sample in samples) {
  dotchart(sample, main=title2, xlab=x.axis)
}

# Histogram of locations of palindromes in original data and uniform scatter
bins <- 35
hist(locations, col=color, breaks=bins, main=title1, xlab=x.axis)

for (sample in samples) {
  hist(sample, breaks=bins, main=title2, xlab=x.axis)
}

# Scatterplot of spacing between consecutive palindromes
title1 <- 'Spacing between Consecutive Palindromes'
title2 <- c(title1,'(Simulated)')
x.axis <- 'Base Pair Location'
y.axis <- 'Distance (Base Pairs) from Previous Palindrome'
y.range <- c(0,5000)
plot(locations[-1], diff(locations), col=color, ylim=y.range, main=title1, xlab=x.axis, ylab=y.axis)

for (sample in samples) {
  plot(sample[-1], diff(sample), ylim=y.range, main=title2, xlab=x.axis, ylab=y.axis)
}

# Histogram of counts of palindromes in non-overlapping regions in original data and uniform scatter
interval.length <- 2500
title1 <- paste('Number of Palindromes in Non-Overlapping Regions of Length', interval.length)
title2 <- c(title1,'(Simulated)')
x.axis <- 'Number of Palindromes'
bins <- seq(0,20,1)
hist(as.vector(table(cut(locations, breaks=seq(0,N,interval.length), include.lowest=TRUE))), breaks=bins, col=color, main=title1, xlab=x.axis)

for (sample in samples) {
  hist(as.vector(table(cut(sample, breaks=seq(0,N,interval.length), include.lowest=TRUE))), breaks=bins, main=title2, xlab=x.axis)
}

Scenario 2: Locations and Spacings

Use graphical methods to examine the spacings between consecutive palindromes and sum of consecutive pairs, triplets, etc, spacings. Compare what you find to what you would expect to find in a random scatter. Also, use graphical methods to compare locations of the palindromes.

# Chi-square Goodness of Fit Test
# Case 1: k(number of sub-intervals)=20
k <- 20
locations.expected <- n/k
tab <- table(cut(locations, breaks=seq(0, N, length.out=k+1), include.lowest=TRUE))
locations.observed <- as.vector(tab)
chi_2 <- sum((locations.observed - locations.expected)^2/locations.expected)
chi2_compare <- qchisq(p=0.95, df=19)
p_value <- pchisq(chi_2, df=19, lower.tail=FALSE)
print(cat('\nWhen conducting chi_square Goodness of fit test comparing locations(divided in 20 sub-intervals) against uniform distribution\n'))

When conducting chi_square Goodness of fit test comparing locations(divided in 20 sub-intervals) against uniform distribution
NULL

print(paste('The value of chi_square statistic is', chi_2))

[1] "The value of chi_square statistic is 17.9189189189189"

print(paste('The p_value is', p_value))

[1] "The p_value is 0.527860332119311"

## Visualization of the Residual
Residuals <- (locations.observed - locations.expected) / sqrt(locations.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 20 sub-intervals)')

# Case 2: k(number of sub-intervals)=30
k <- 30
locations.expected <- n/k
tab <- table(cut(locations, breaks=seq(0, N, length.out=k+1), include.lowest=TRUE))
locations.observed <- as.vector(tab)
chi_2 <- sum((locations.observed - locations.expected)^2/locations.expected)
chi2_compare <- qchisq(p=0.95, df=29)
p_value <- pchisq(chi_2, df=29, lower.tail=FALSE)
print(cat('\nWhen conducting chi_square Goodness of fit test comparing locations(divided in 30 sub-intervals) against uniform distribution\n'))

When conducting chi_square Goodness of fit test comparing locations(divided in 30 sub-intervals) against uniform distribution
NULL

print(paste('The value of chi_square statistic is', chi_2))

[1] "The value of chi_square statistic is 40.6891891891892"

print(paste('The p_value is', p_value))

[1] "The p_value is 0.0732835870345071"

## Visualization of the Residual
Residuals <- (locations.observed - locations.expected) / sqrt(locations.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 30 sub-intervals)')

# Case 3: k(number of sub-intervals)=60
k <- 60
locations.expected <- n/k
tab <- table(cut(locations, breaks=seq(0, N, length.out=k+1), include.lowest=TRUE))
locations.observed <- as.vector(tab)
chi_2 <- sum((locations.observed - locations.expected)^2/locations.expected)
chi2_compare <- qchisq(p=0.95, df=59)
p_value <- pchisq(chi_2, df=59, lower.tail=FALSE)
print(cat('\nWhen conducting chi_square Goodness of fit test comparing locations(divided in 60 sub-intervals) against uniform distribution\n'))

When conducting chi_square Goodness of fit test comparing locations(divided in 60 sub-intervals) against uniform distribution
NULL

print(paste('The value of chi_square statistic is', chi_2))

[1] "The value of chi_square statistic is 79"

print(paste('The p_value is', p_value))

[1] "The p_value is 0.0421403871302519"

## Visualization of the Residual
Residuals <- (locations.observed - locations.expected) / sqrt(locations.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 60 sub-intervals)')

# Histogram of locations of palindromes in original data and uniform scatter
sample <- runif(n, min=0, max=N)
title <- 'Locations of Palindromes (Original vs. Simulated)'
x.axis <- 'Base Pair'
bins <- 35
hist(locations, breaks=bins, probability=TRUE, col=rgb(1,0,0,0.5), main=title, xlab=x.axis)
lines(density(locations, adjust=2), col=2)

hist(sample, breaks=bins, probability=TRUE, col=rgb(0,0,1,0.5), add=TRUE)
lines(density(sample, adjust=2), col=4)

legend('topright', legend=c('Original', 'Uniform'), lty=c(1,1), col=c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Single Palindrome Spacing
locations.sorted = sort(locations, decreasing = FALSE)
distance.single <- abs(locations.sorted[-1]-locations.sorted[-length(locations.sorted)])
# Histogram of spacings of palindromes in original data and exponential distribution
hist(distance.single, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Consecutive palindrome Spacings Distribution Comparison", xlab = "Distance between Consecutive Palindromes", ylim = c(0,0.001))
lines(density(distance.single, adjust = 2), col = rgb(1,0,0,0.5))

Expo <- rexp(n-1, rate = 1/mean(distance.single))
hist(Expo, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Expo, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 4200, y = 0.0009, legend = c("Sample", "Exponential"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-square Goodness of Fit Test
# Case 1: Divided in 7 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.single, decreasing = FALSE)
lambda <- 1/mean(distance.single)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(0,0.05, 0.1, 0.3, 0.5, 0.7, 0.9,1)))
spacings.expected <- (n-1)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.single, breaks=spacings.intervals, include.lowest=TRUE)))
contingency_7 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_7

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=5)
p_value <- pchisq(chi_2, df=5, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 7 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 7 sub-intervals is 6.42447875347949e-05"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 7 bins)')

# Case 2: Divided in 10 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.single, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.1))))
spacings.expected <- (n-1)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.single, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_10 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_10

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=8)
p_value <- pchisq(chi_2, df=8, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 10 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 10 sub-intervals is 0.000218981752365422"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 10 bins)')

# Case 3: Divided in 20 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.single, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.05))))
spacings.expected <- (n-1)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.single, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_20 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_20

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=18)
p_value <- pchisq(chi_2, df=18, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 20 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 20 sub-intervals is 0.0117020696134169"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Palindrome locations', main='Plot of Standardized Residual for Locations (divided in 20 bins)')

# Consecutive Pairs
locations.sorted <-  sort(locations, decreasing = FALSE)
locations.pairs <- locations.sorted[-length(locations.sorted)]
distance.pairs <- abs(locations.sorted[-1][-1]-locations.pairs[-length(locations.pairs)])
# Histogram of spacings of palindromes in original data and exponential distribution
hist(distance.pairs, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Consecutive Pairs Spacings Distribution Comparison", xlab = "Distance between Consecutive Pairs of Palindromes Locations", ylim = c(0,0.001))
lines(density(distance.pairs, adjust = 2), col = rgb(1,0,0,0.5))

Expo <- rexp(n-2, rate = 1/mean(distance.pairs))
hist(Expo, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Expo, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 4200, y = 0.0005, legend = c("Sample", "Exponential"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-square Goodness of Fit Test
# Case 1: Divided in 7 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.pairs, decreasing = FALSE)
lambda <- 1/mean(distance.pairs)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(0,0.05, 0.1, 0.3, 0.5, 0.7, 0.9,1)))
spacings.expected <- (n-2)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_7 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_7

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=5)
p_value <- pchisq(chi_2, df=5, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 7 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 7 sub-intervals is 4.0192806354887e-05"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Spacings between Palindrome Pairs', main='Plot of Standardized Residual for Locations (divided in 7 bins)')

# Case 2: Divided in 10 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.pairs, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.1))))
spacings.expected <- (n-2)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_10 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_10

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=8)
p_value <- pchisq(chi_2, df=8, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 10 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 10 sub-intervals is 1.68753795541538e-05"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Spacings between Palindrome Pairs', main='Plot of Standardized Residual for Locations (divided in 10 bins)')

# Case 3: Divided in 20 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.pairs, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.05))))
spacings.expected <- (n-2)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_20 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_20

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p=0.95, df=18)
p_value <- pchisq(chi_2, df=18, lower.tail=FALSE)
print(paste('The p_value when the distance is splited into 20 sub-intervals is', p_value))

[1] "The p_value when the distance is splited into 20 sub-intervals is 0.000159225622786541"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Spacings between Palindrome Pairs", main = "Plot of Standardized Residual for Locations (divided in 20 bins)")

# Consecutive Triplets
locations.sorted <-  sort(locations, decreasing = FALSE)
locations.triplets <- locations.sorted[-length(locations.sorted)]
locations.triplets <- locations.triplets[-length(locations.triplets)]
distance.triplets <- abs(locations.sorted[-1][-1][-1]-locations.triplets[-length(locations.triplets)])
# Histogram of spacings of palindromes in original data and exponential distribution
hist(distance.triplets, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Consecutive Triplets Spacings Distribution Comparison", xlab = "Distance between Consecutive Palindromes Triplets", ylim = c(0,0.0004))
lines(density(distance.triplets, adjust = 2), col = rgb(1,0,0,0.5))

Gam <- rgamma(n-2, 2, rate = 1/mean(distance.pairs))
hist(Gam, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Gam, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 5000, y = 0.0003, legend = c("Sample", "Gamma"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-square Goodness of Fit Test (Need to be changed)
# Case 1: Divided in 7 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.triplets, decreasing = FALSE)
lambda <- 2/mean(distance.pairs)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(0,0.05, 0.1, 0.3, 0.5, 0.7, 0.9,1)))
spacings.expected <- (n-3)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_7 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_7

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p = 0.95, df = 5)
p_value <- pchisq(chi_2, df = 5, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 7 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 7 sub-intervals is 9.07421505828639e-47"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome locations", main = "Plot of Standardized Residual for Locations (divided in 7 bins)")

# Case 2: Divided in 10 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.pairs, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.1))))
spacings.expected <- (n-3)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_10 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_10

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p = 0.95, df = 8)
p_value <- pchisq(chi_2, df = 8, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 10 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 10 sub-intervals is 8.252676182357e-37"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome locations", main = "Plot of Standardized Residual for Locations (divided in 10 bins)")

# Case 3: Divided in 20 intervals
# Construct expected number of intervals
spacings.observed <- sort(distance.pairs, decreasing = FALSE)
spacings.intervals <- as.numeric(quantile(spacings.observed, probs = c(seq(0,1, by = 0.05))))
spacings.expected <- (n-3)*(exp(-lambda*spacings.intervals[-length(spacings.intervals)])-exp(-lambda*spacings.intervals[-1]))
spacings.expected[length(spacings.expected)] <- n-sum(spacings.expected[1:length(spacings.expected)-1])
spacings.observed <- as.numeric(table(cut(distance.pairs, breaks = spacings.intervals, include.lowest = TRUE)))
contingency_20 <- data.frame(spacings.intervals[-1],spacings.observed,spacings.expected)
contingency_20

chi_2 <- sum((spacings.observed - spacings.expected)^2/spacings.expected)
chi2_compare <- qchisq(p = 0.95, df = 18)
p_value <- pchisq(chi_2, df = 18, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 20 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 20 sub-intervals is 1.34421275018616e-39"

## Visualization of the Residual
Residuals <- (spacings.observed - spacings.expected) / sqrt(spacings.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome locations", main = "Plot of Standardized Residual for Locations (divided in 20 bins)")

Scenario 3: Counts

Use graphical methods and more formal statistical tests to examine the counts of palindromes in various regions of the DNA. Split the DNA into nonoverlapping regions of equal length to compare the number of palindomres in an interval to the number of that you would expect from uniform random scatter. The counts for shorter regions will be more variable than those for logner regions. Also, consider classifying the regions according to the number of counts.

regionsplit <- function(n.region, gene, site){
  count.int <- table(cut(site, breaks = seq(1, length(gene), length.out=n.region+1), include.lowest=TRUE))
  count.vector <- as.vector(count.int)
  count.tab <- table(count.vector)
  return (count.tab)
}

# Case 1: divided by 40 intervals
n.region <- 40
gene <- seq(1,N)
observed <- as.numeric(regionsplit(n.region, gene, locations))
interval <- as.numeric(names(regionsplit(n.region, gene, locations)))
lambda <- n/n.region
# Histogram of counts of palindromes in original data and poisson distribution
counts.hist <- regionsplit(n.region, gene, locations)
hist(counts.hist, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Counts of Palindrome Distribution Comparison (40 sub-intervals)", xlab = "Counts found per interval", ylim = c(0,0.6))
lines(density(counts.hist, adjust = 2), col = rgb(1,0,0,0.5))

Pois <- rpois(n, lambda)
hist(Pois, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Pois, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 7, y = 0.5, legend = c("Sample", "Poisson"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-sqr Goodness of Fit test
expected <- n.region*exp(-lambda)* lambda**(interval)/factorial(interval)
for (i in c(0:2)){
  expect <- n.region*exp(-lambda)* lambda**(i)/factorial(i)
  expected[1] <- expected[1]+ expect
}
expect <- n.region*exp(-lambda)* lambda**(11)/factorial(11)
expected[8] <- expected[8]+ expect
expected[10] <- 0
for (i in c(1:12)){
  expect <- exp(-lambda)* lambda**(i)/factorial(i)
  expected[10] <- expected[10]+ expect
}
expected[10] <- (1-expected[10])*n.region
counts.expected <- c()
counts.interval <- c()
counts.observed <- c()
# Group bins
counts.expected[1] <- sum(expected[1:2])
counts.expected[2] <- sum(expected[3:4])
counts.expected[3] <- sum(expected[5])
counts.expected[4] <- sum(expected[6:7])
counts.expected[5] <- sum(expected[8:10])
counts.observed[1] <- sum(observed[1:2])
counts.observed[2] <- sum(observed[3:4])
counts.observed[3] <- sum(observed[5])
counts.observed[4] <- sum(observed[6:7])
counts.observed[5] <- sum(observed[8:10])
counts.interval[1] <- interval[2]
counts.interval[2] <- interval[4]
counts.interval[3] <- interval[5]
counts.interval[4] <- interval[7]
counts.interval[5] <- interval[7]+1
counts.table40 <- data.frame(counts.interval,counts.observed,counts.expected)
counts.table40

# Chi-square statistic
chi_2 <- sum((counts.observed - counts.expected)^2/counts.expected)
chi2_compare <- qchisq(p = 0.95, df = 3)
p_value <- pchisq(chi_2, df = 3, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 40 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 40 sub-intervals is 0.365205436984496"

## Visualization of the Residual
Residuals <- (counts.observed - counts.expected) / sqrt(counts.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome counts", main = "Plot of Standardized Residual for Counts (divided in 40 sub-intervals)")

# Case 2: divided by 60 intervals
n.region <- 60
gene <- seq(1,N)
observed <- as.numeric(regionsplit(n.region, gene, locations))
interval <- as.numeric(names(regionsplit(n.region, gene, locations)))
lambda <- n/n.region
# Histogram of counts of palindromes in original data and poisson distribution
counts.hist <- regionsplit(n.region, gene, locations)
hist(counts.hist, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Counts of Palindrome Distribution Comparison (60 sub-intervals)", xlab = "Counts found per interval", ylim = c(0,0.6))
lines(density(counts.hist, adjust = 2), col = rgb(1,0,0,0.5))

Pois <- rpois(n, lambda)
hist(Pois, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Pois, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 7, y = 0.5, legend = c("Sample", "Poisson"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-sqr Goodness of Fit test
expected <- n.region*exp(-lambda)* lambda**(interval)/factorial(interval)
expect <- n.region*exp(-lambda)* lambda**(0)/factorial(0)
expected[1] <- expected[1]+ expect
for (i in c(10:11)){
  expect <- n.region*exp(-lambda)* lambda**(i)/factorial(i)
  expected[9] <- expected[9]+ expect
}
expected[11] <- 0
for (i in c(1:12)){
  expect <- exp(-lambda)* lambda**(i)/factorial(i)
  expected[11] <- expected[11]+ expect
}
expected[11] <- (1-expected[11])*n.region
counts.expected <- c()
counts.interval <- c()
counts.observed <- c()
# Group bins
counts.expected[1] <- sum(expected[1:2])
counts.expected[2:6] <- expected[3:7]
counts.expected[7] <- sum(expected[8:11])
counts.observed[1] <- sum(observed[1:2])
counts.observed[2:6] <- observed[3:7]
counts.observed[7] <- sum(observed[8:11])
counts.interval[1] <- interval[2]
counts.interval[2:6] <- interval[3:7]
counts.interval[7] <- interval[7]+1
counts.table60 <- data.frame(counts.interval,counts.observed,counts.expected)
counts.table60

# Chi-square statistic
chi_2 <- sum((counts.observed - counts.expected)^2/counts.expected)
chi2_compare <- qchisq(p = 0.95, df = 5)
p_value <- pchisq(chi_2, df = 5, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 60 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 60 sub-intervals is 0.406748395085584"

## Visualization of the Residual
Residuals <- (counts.observed - counts.expected) / sqrt(counts.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome counts", main = "Plot of Standardized Residual for Counts (divided in 60 sub-intervals)")

# Case 3: divided by 80 intervals
n.region <- 80
gene <- seq(1,N)
observed <- as.numeric(regionsplit(n.region, gene, locations))
interval <- as.numeric(names(regionsplit(n.region, gene, locations)))
lambda <- n/n.region
# Histogram of counts of palindromes in original data and poisson distribution
counts.hist <- regionsplit(n.region, gene, locations)
hist(counts.hist, breaks= 15, col = rgb(1,0,0,0.5), probability = TRUE, main = "Counts of Palindrome Distribution Comparison (80 sub-intervals)", xlab = "Counts found per interval", ylim = c(0,0.6))
lines(density(counts.hist, adjust = 2), col = rgb(1,0,0,0.5))

Pois <- rpois(n, lambda)
hist(Pois, breaks = 15, col = rgb(0,0,1,0.5), probability = TRUE, add = TRUE)
lines(density(Pois, adjust = 2), col = rgb(0,0,1,0.5))

legend(x = 7, y = 0.5, legend = c("Sample", "Poisson"), lty = c(1,1), col = c(rgb(1,0,0,0.5), rgb(0,0,1,0.5)))

# Chi-sqr Goodness of Fit test
expected <- n.region*exp(-lambda)* lambda**(interval)/factorial(interval)
for (i in c(8:9)){
  expect <- n.region*exp(-lambda)* lambda**(i)/factorial(i)
  expected[9] <- expected[9]+ expect
}
expected[10] <- 0
for (i in c(1:10)){
  expect <- exp(-lambda)* lambda**(i)/factorial(i)
  expected[10] <- expected[10]+ expect
}
expected[10] <- (1-expected[10])*n.region
counts.expected <- c()
counts.interval <- c()
counts.observed <- c()
# Group bins
counts.expected[1] <- sum(expected[1:2])
counts.expected[2:6] <- expected[3:7]
counts.expected[7] <- sum(expected[8:10])
counts.observed[1] <- sum(observed[1:2])
counts.observed[2:6] <- observed[3:7]
counts.observed[7] <- sum(observed[8:10])
counts.interval[1] <- interval[2]
counts.interval[2:6] <- interval[3:7]
counts.interval[7] <- interval[7]+1
counts.table60 <- data.frame(counts.interval,counts.observed,counts.expected)
counts.table60

# Chi-square statistic
chi_2 <- sum((counts.observed - counts.expected)^2/counts.expected)
chi2_compare <- qchisq(p = 0.95, df = 5)
p_value <- pchisq(chi_2, df = 5, lower.tail = FALSE)
print(paste("The p_value when the distance is splited into 60 sub-intervals is", p_value))

[1] "The p_value when the distance is splited into 60 sub-intervals is 0.868215562053925"

## Visualization of the Residual
Residuals <- (counts.observed - counts.expected) / sqrt(counts.expected)
plot(Residuals, type = 'h', ylab = "Standardized Residuals", xlab = "Palindrome counts", main = "Plot of Standardized Residual for Counts (divided in 80 sub-intervals)")

Scenario 4: The Biggest Cluster

Does the interval with the greatest number of palindromes indicate a potential origin of replication? Be careful in making your intervals, for any small, but significant deviations from random scatter, such as a tight cluster of a few palindromes, could easily go undetected if the regions examined are too large. Also, if the regions are too small, a cluster of palindromes may be split between adjacent interavls and not appear as a high-count interval.

final <- array(dim=c(500,1))
interval_length <- array(dim=c(500,1))
lamda <- array(dim=c(500,1))
for (k in 20:100){
  tab <- table(cut(locations, breaks=seq(0, N, length.out=k+1), include.lowest=TRUE))
  head(tab,10)
  tab<-as.vector(tab)
  lamda[k,] <-sum(tab)/k
  threshold <-max(tab)
  result <- 0
  interval_length[k,] <- N/k
  for (i in 0:(threshold-1)){
    result <- result+((lamda[k]^i)*exp(-lamda[k])/factorial(i))
  }
  final[k,] <- 1-result^k
}
result <- data.frame(lamda,interval_length,final)
# Display Table containing the probability of a Poisson Distribution having e greatest number of hits at least k for each sub-interval divisions
result[c(40,60,80),]

Additional Scenario: HIV and Age

TODO Description

# Clean out 'unknown' data and convert factor to numerical values
health <- transform(health, age_yrs=as.numeric(age_yrs),
                            hiv=as.character(hiv))
health.ind <- which(health$hiv != 'unknown')
health <- health[health.ind,]
# Total number of people that have hiv
population=nrow(health)
pop_hiv <- nrow(health[which(health$hiv=='positive'),])
# Split the age into four groups
# 0-20
age_first <- health$age_yrs[which(health$age_yrs<21)]
age_proportion_first <- length(age_first)/population
hiv_proportion_first<- nrow(health[which((health$hiv== 'positive') & (health$age_yrs <21)),])/pop_hiv
# 21-40
age_second<-health$age_yrs[which(health$age_yrs>20 & health['age_yrs']<41)]
age_proportion_second <- length(age_second)/population
hiv_proportion_second<-nrow(health[which(health$age_yrs>20 &health$age_yrs<41 & health$hiv=='positive'),])/pop_hiv
# 41-60
age_third<-health$age_yrs[which(health['age_yrs']>40 & health['age_yrs']<61)]
age_proportion_third <- length(age_third)/population
hiv_proportion_third<-nrow(health[which(health$age_yrs>40 & health$age_yrs<61 &health$hiv=='positive'),])/pop_hiv
# 61+
age_last<-health$age_yrs[which(health['age_yrs']>60)]
age_proportion_last <- length(age_last)/population
hiv_proportion_last<-nrow(health[which(health$age_yrs>60 & health$hiv=='positive'),])/pop_hiv
# Expected Data
population_dist <-c(age_proportion_first,age_proportion_second,age_proportion_third,age_proportion_last)
# Observed Data
hiv_dist<-c(hiv_proportion_first,hiv_proportion_second,hiv_proportion_third,hiv_proportion_last)
age_dist <- c("0-20", "21-40", "41-60", "61+")
data.frame(age_dist,population_dist,hiv_dist)

# Goodness-fittest
chi_2 <- sum((hiv_dist - population_dist)^2/population_dist)
chi2_compare <- qchisq(p=0.95, df=3)
p_value <- pchisq(chi_2, df=3, lower.tail=FALSE)
print(paste('The p_value of Goodness of Fit Test is',p_value))

[1] "The p_value of Goodness of Fit Test is 0.919175049520472"

#Visualization
Residuals <- (hiv_dist - population_dist) / sqrt(population_dist)
plot(Residuals, type='h', ylab='Standardized Residuals', xlab='Proportion of Positive HIV', main='Plot of Standardized Residual for Age and HIV Positive')

Null Hypothesis: The proportion of age in the population is unrelated with the proportion of people having hiv.(Age is not an influencing factor for HIV testing positive) Since p-value of this chi-square goodness of fit test is close to 1, we see that deviations as large as ours (or larger) are very likely. In addition, having values of the standardized residual less than 3 suggests that it is a good fit of the age distribution to estimate the people testing positive on hiv. Hence, we reject the null hypothesis and conclude that the distribution of proportion of age matches with the the distribution of people testing positive on HIV.